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ABSTRACT 

q 

■rj" , A new technique for reliably identifying point sources in millimeter /sub-millimeter 

!^ \ wavelength maps is presented. This method accounts for the frequency dependence 

of noise in the Fourier domain as well as non-uniformities in the coverage of a field. 
This optimal filter is an improvement over commonly-used matched filters that ignore 

►v> ■ coverage gradients. Treating noise variations in the Fourier domain as well as map 

H ' space is traditionally viewed as a computationally intensive problem. We show that 

■ ■ ■ the penalty incurred in terms of computing time is quite small due to casting many of 

the calculations in terms of FFTs and exploiting the absence of sharp features in the 
noise spectra of observations. Practical aspects of implementing the optimal filter are 
presented in the context of data from the AzTEC bolometer camera. The advantages 
of using the new filter over the standard matched filter are also addressed in terms of a 
typical AzTEC map. 

Subject headings: Data Analysis and Techniques, Astrophysical Data 



Introduction 



The discovery and study of sub-millimeter galaxies, or SMGs, has become a key enterprise 
wit hin millimeter /su b- millimeter (mm/sub-mm) astronomy over the past decade or so (see review 
by iBlain et alj |2002| ). With bolometric luminosities ^ 5 x 10^^ Lq and star formation rates > 
lOOM0yr~^ , SMGs represent some of the most luminous galaxies in the early Universe. The 
strong negative /c-correction at A > 500 /xm means that a galaxy of a given luminosity will be 
equally detectable in a flux-limited survey from 1 ^ -z ^ 10. A deep, wide-area survey at mm/sub- 
mm wavelengths is thus a sensitive probe of starburst galaxies from the epoch at which they first 
turn on through the peak of star formation activity in the Universe at z ~ 1 — 2. Several papers 
have reported on th e source counts of SMGs detected at A = 250 — 2000 Mm (e.g. lCoppin et al.ll2006l : 
Bertoldi et aUboO?! : Iweifi et al-lbood : Ivieira et allboiol : IScott et ahlboid . and refere nces therein) 



which provide strong cons t raints on the evo l utioii ary history of massive galaxies (e.g. iBaugh et al 



20051 : IValiante et al.ll2009l : iBethermin et al.ll2012l ). SMG surveys from single-dish telescopes also 
provide a catalog of interesting targets to follow-up with higher angular resolution imaging of the 
dust and molecular gas in these objects, which inform on the physical processes that trigger and 
maintain the starbursts in these galaxies. 

An important criterion for the success of such surveys is the design of a data analysis scheme 
that can reliably identify SMG candidates on sky maps. It is especially important to minimize the 
rate of false detections because they can result in a waste of valuable time and resources in follow- 
up observations. But reliably identifying astronomical signal in mm/sub-mm maps is a difficult 
task because these maps are inherently low in signal/noise due to atmospheric contamination, 
instrument noise, and the high confusion limit of surveys from typical mm/sub-mm telescopes 
(~10m diameter). In the wide field surveys designed to detect SMGs, they usually appear as 
point sources much smaller than the angular resolution of the telescopes used. A simple appro ach 
to identify point sources that is used by standard data analysis packages (see IStetsonI 119871 . for 
example) is: 1) smooth the signal map by convolving it with the point spread function (PSF); 
2) assume that the errors associated with pixels are uncorrelated and propagate them through 
the convolution; and 3) select the peaks in the smoothed map based on the per-pixel signal/noise 
determined from steps 1) and 2). This reduces high-frequency signal variations between adjacent 
pixels and thus increases the signal/noise for point-source detection. 

However, there are two general problems with real mm/sub-mm maps that renders the above 
procedure inadequate for reliably identifying point sources. First, the noise in these maps is gen- 
erally not white; it is more pronounced on larger scales or at low (spatial) frequencies. This effect, 
which gives rise to pixel-pixel correlations, is due mainly to 1// drifts in atmospheric and/or in- 
strumental conditions. In mm/sub-mm astronoi ny, this problem is usually treated with a "matched 
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filter " implemented in the Fourier domain (e.g. iTegmark &: de 01iveira-Costalll998l : iBarreiro et al. 
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However, this 

filtering technique is optimal only in the case of uniform coverage, which brings up the second prob- 
lem with real mm/sub-mm maps: the coverage of a field is non-uniform in general. Because these 



observations are usually carried out by scanning an array of detectors across the field of interest, 
variations in atmospheric conditions or detector noise during the scan are often responsible for this 
non-uniformity. With commonly used schemes such as raster- or Lissajous-scanning, coverage also 
tends to decrease smoothly from the center to the edges of the map. 

The correct way to deal with both these problems involves the construction and inversion of 
a pixel-pixel noise covariance matrix. This path can be extremely challenging computationally for 
maps containing ~ 100, 000 pixels. Therefore, a common strategy is to pick out a region of the map 
that has essentially uniform coverage and then apply the standard matched filter to it. Even when 
a near-uniform region exists, ignoring small coverage variations within it can have a noticeable 
effect, as we will show below. Furthermore, despite gradients, the coverage may be deep enough 
near the edges of the map to identify bright sources with high significance, and these "border" 
regions can often cover a significant area compared to the near-uniform region. In some cases, 
significant coverage gradients are unavoidable due to difficult observing conditio ns or an observing 



strategy where a mosaic of small maps are stitched together (for example, see iBorys et al. 
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Austermann et al.l l2010l ) . A truly optimal analysis would help improve the 



reliability with which point-source candidates are identified in all regions of a map despite coverage 
gradients. 

In this paper, we present an algorithm that addresses non-uniform coverage as well as low- 
frequency noise in a computationally efficient manner. This technique is designed to be optimal in 
the regime where blending of resolved sources is negligible. The work presented here is an extension 
of the standard data reduction tool s used for point-sou rce extraction in maps taken with the 1 . 1 mm 
bolometer array c amera, AzTEC (IWilson et al.ll2008l) . While this builds on the existing AzTEC 



reduction pipeline (jScott et al. 



2008 
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20081 ). the principles are generic to any scan-data 



from mm/sub-mm arrays. In section [2] we present the formal methodology that will be used to 
identify point sources and to estimate their brightness. Essentially, this method involves a least- 
squares fit of the PSF to every point on the map. In section [Sj we demonstrate that, for the case of 
uniform coverage, the fit paradigm converges to the familiar and computationally efficient matched 
filter routinely used in point-source searches. In section HI we develop an implementation of the 
pixel-by-pixel fit that can efficiently handle the case of non-uniform coverage. Then, in section \5\ 
we discuss the practicalities of applying this technique to AzTEC data. Finally, in section [6l we 
conclude by discussing the performance and scope of the optimal filter introduced here. 



2. The Method 



We leverage our knowledge of the shape of isolated point sources in our maps to construct an 
optimal or wiener filter using the point spread function (PSF) as the template. Applying this filter 
to the map is formally equivalent to centering the PSF on each pixel of the map and fitting for 
the best-fit amplitude (IStetsonlll987l ). In our case, we store the best fit amplitude at each pixel in 
a separate map which we will refer to as the filtered signal map. In addition, we construct a new 



map of the same dimensions as the signal map which contains 1/error^ estimates of the best fit 
amphtudes. This map will be referred to as the filtered coverage map. These maps, augmented by 
noise realizations of the field, are the primary inputs to the identification of point sources. 



2.1. AzTEC data as an example 



AzTEC is a 144-element semiconductor- type bolometer array that imaged the sky at 1.1-mm 
wavelength over two successful observing cam paigns: one using the 15-m James Clark Maxwell 
Telescope (JCMT) in Hawaii from 2005-2006 (jwilson et aljboosl ). and the other using the 10-m 
Ataca ma Sub-millimeter Telescope Experiment (ASTE) in the Atacama desert of Chile from 2007- 
2008 (JEzawa et al.ll2008l ). Over these two observing runs, close to a hundred fields were mapped, 
and each map comprises 100,000 pixels or more. Therefore, computational speed was an important 
consideration when developing this algorithm. Computational speed will be even more important 
for newer and upcoming mm/sub-mm observatories like the Large Millimeter Telescope and CCAT 
that can image the sky at even faster rates with higher angular resolution. 






157.7 157.6 

RA [dec] 



57.5 157.4 



mJy 



157.7 157.6 

RA [deq] 



(a) 



(b) 



Fig. 1. — (a) 1.1mm- wavelength map of a region surrounding the quasar SDSSJ1030-I-0524, ob- 
tained with AzTEC- ASTE. The outer regions are saturated high (black) or low (white) on an 
astronomically interesting scale (itSmJy). Astronomical features are visible in the central un- 
saturated (grey scale) region, mainly as PSF-sized dark patches ("blobs"), or point sources, (b) 
Coverage map of this field. It is filled with the values l/e^(xp), where e(xp) is the uncertainty of 
the signal at the p'th pixel of (a). Coverage contours are labeled relative to the maximally covered 
(dark) central regions. 



To illustrate the starting point for the filtering process, we present in Fig. [T] an unfiltered 
AzTEC signal map and the corresponding coverage map. In the AzTEC data analysis, coverage 
maps have units of 1/Jy^ and are filled with estimates of 1/error^ at each pixel. In Fig. [1] and 
throughout the rest of this paper, we will use maps and plots related the AzTEC-ASTE data on 
a region centered around the high-redshift quasar SDSSJ1030+0524 to illustrate our points. This 
data set consists of 45 observations of the field, each lasting about forty minutes, carried out over 



Novem ber 2008. Sc i entific re sults from this AzTEC-ASTE field are published in iHumphrey et al. 



(|201ll ) and lZeballosI (jin prep. I ). These data were taken by scanning the AzTEC array in a Lissajous 
pattern, and results in a map with deepest coverage in the center, and decreasing coverage toward 
the edges of the map. 



2.2. Assumptions 

Two assumptions, both of which are generally applicable to most observations, are critical for 
the efficient implementation of the process described below. They are: 

1. The effective PSF does not change with location on the map. We represent the PSF by the two- 
dimensional function /(x). The value of this function at the z'th pixel when the PSF is centered 
on the p'th pixel will be denoted by /(xi — Xp), where Xi and Xp are 2-d vectors that specify the 
locations of the i'th and p'th pixels, respectively. We note that while this assumption is explicitly 
not valid for some instruments (e.g., the Chandra X-ray telescope), this is a reasonable assumption 
for many imaging instruments. 

2. Noise-induced correlations between two pixels, hereafter denoted Corr(xk,xi), depend only on 
the distance |xk — xi|, between the pixels. That is, neither the location of the two pixels on the map 
nor their relative orientation determines the correlation in their noise. In such cases, Corr(xj5;,xi) 
may be conveniently expressed as the inverse Fourier transform (IFT) of the power spectral density 
(PSD) of noise in the map. If we denote the PSD by y^(ka), where ka is a 2-d vector in the Fourier 
domain of the map, 



N, 



pixel 



Corr(xk,xi) = ^ l^^(ka) exppvrjka • (xk - xi)], (1) 

a 

where A'^pixei is the total number of pixels in the map, and the sum is over all A'^pixei vectors ka- 
(Appendix A contains a brief description of the discrete Fourier transform conventions and their 
corollaries used here.) 

In Eq. m y^(ka) is a normalized form of the PSD that satisfies Zl^"''"'' ^^(ka) = 1, so that 
the diagonal elements of the correlation matrix evaluate to 1. Note that for a completely flat PSD 
(y^(ka) = cnst) the correlation matrix is diagonal. In most observations the noise PSD varies 
smoothly with k and has broad features rather than narrowly peaked ones (see Fig. [3l for example). 
In such cases the correlation matrix is band-diagonal (i.e., Corr(xk,xi) ~ when |xk — xi| is large). 



This is a property that we will exploit below in our pursuit of computational efficiency. 



2.3. An optimal filter via a generalized least-squares fit 

We implement the optimal filter by minimizing the quantity 

N ■ 1 

^ "pixel 

xl= J2 ^^(^^) ~ Sp/(^k - Xp)]Wfc/[d(xi) - Sp/(xi - Xp)], (2) 

k,l=0 

for each pixel, p, in the map where (i(xk) and (i(xi) denote the value of the unfiltered signal map 
at the fc'th and /'th pixels, Sp is the amplitude of the fitted PSF, and W is an A'^pixei x A'pixci weight 
matrix whose calculation is described below. The summations in Eq. [2] are over all pixels of the 
map. Eq. [2]is minimized when 



^W^feid(xk)/(xi-Xp) 



kl 



^ Wkifi^k - Xp)/(xi - Xp) 

k,l 



(3) 



For least-squares fits, the minimum variance choice of the weight matrix is W = C~^ where 
C is the pixel-pixel noise covariance matrix (i.e. the map's covariance in the absence of signal). In 
this case, the Xp calculated in Eq. [2] is drawn from a true x^ distribution with A'pixcis — 1 degrees 
of freedom. The elements of C are 



Cm = 6(xk)Corr(xk,xi)6(xi) = e(xk) | J^ y2(]^^)^2.,k..(x,-x,) J ^(^^) 



(4) 



where e(xk) and e(xi) are the standard deviations of noise in the k-th and l-th pixels respectively. 
Note that the diagonal elements, Ckk evaluate to e^(xk), as they should. 

The elements of the weight matrix can then be calculated as 

Using Eq. IA5I of Appendix A, it is straightforward to verify that J2i ^kiCim = ^km- Substituting 
the above form of Wm into Eq. [3l we may express the numerator and denominator of Eq. [3] as 



AT - 1 ^^ / (Xk - Xp) ^ ,^ g.^,.a-,x.-xi; ^ ^^^^^ 



n _ _J_ sr /i^k-xpj ^ e-'^---^-'— ^ /(xi-xpj 



In practice, we form separate "maps" corresponding to A^ and D because the latter is useful by 
itself (see below). In doing so, the nested sums of Eq. [6] and Eq. [7] will be recognized as Fourier 
transforms, allowing us to benefit from the efficiency of the EFT algorithm. The element-by-element 
ratio of these two, Sp = Np/Dp, forms the filtered signal map. 



2.4. Error propagation 

The error in our estimate of Sp, which we denote Up, is the rms deviation of the best fit PSE 
amplitude, Sp, from the actual amplitude, a^, of a point source that is centered on the p'th pixel, 
i.e. rip = {\sp — flpP). To estimate rip, we consider the influence of a point source located at the p'th 
map pixel on a neighboring pixel i. 

d(xi) = ap/(xi - Xp) + \fcli. (8) 

As mentioned above, we recognize that the possibility of other resolved sources being quite close to 
pixel i (source blending) is ignored in Eq. [8]and must be treated iteratively. Using Eq. [8]for d(y^\) 
and Eq. [3] for Sp, our estimate for n^ is 

X] X] ^fe«Wf„CH/(xi - Xp)/(Xu - Xp) 

^ VFfci/(xk - Xp)/(xi - Xp 



2 /I i2\ '«/ '-'" icx\ 

n„ = (|sp - Opl ) = -^ , (9) 






and since W = C , the above expression simplifies to 

1 



lY^ Wkifi^k - Xp)/(xi - Xp) 

k,l 



/Dp 



(10) 



Thus, when the optimal W is used, no extra steps are needed to evaluate rip, as Dp already exists 
from generating the filtered signal map. 



2.5. Source significance and goodness of fit 

The filtered signal-to-noise map, filled with values Sp/rip, can be used to both identify and 
locate sources of high significance. Rewriting Eq. [2] in a more convenient form using Eq. [3] and 
Eq. [10] gives 

4 = Y1 Wkid{^^)d{^i) - {sp/ripf. (11) 

k,l 

Because the first term in Eq. [11] is common to all pixels p, the second term is a direct indication of 
the goodness of fit at the p'th pixel. Thus, the locations of peaks in the signal-to-noise map, Sp/rip, 
provide the best estimates of point-source positions on the map. 



3. The case of uniform coverage 

It is instructive to first consider the case of uniform coverage. In this case, e(xk) = e(xi) = cnst 
in Eqs. [6] and [7] and so can be moved outside the summations. The sums over k and / are then 
recognized as Fourier transforms (FTs) of the signal map ((i(x)) and the PSF (/(x)), leading to 
the result that 

/*(ka)d(k. 



E 



^2(k,) 



exppTTJka • Xp 



E 



l/(ka 



(12) 



v^{K 



where /(ka) and d(ka) are the 2-d FTs of /(xj) and (i(xi) respectively. The above expression is the 
familiar "matched filter" that is commonly used in mm/sub - mm astronomy for i dentifying point 
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and is essentially a band-pass filter in the sense that the 
data {d) are low-pass filtered by the PSF / and high-pass filtered by 1/V'^ as the PSD usually 
displays a l//-type trend with spatial frequency. The denominator of Eq. [12] simply provides the 
correct normalization. What Eq. [12] shows is that the point-source finding technique of fitting a 
PSF to each map pixel converges to the conventional matched filter in the case of uniform coverage. 

Another important feature of Eq.[T2]is that, in terms of the FTs /(ka) and (i(ka) the numerator 
and the denominator are each reduced to a single sum of A'pixei terms and the total calculation 
scales as iVpjxei log2 -^pixei rather than as N'^^^^i- As for having to evaluate Sp at each pixel, we note 
that (1) the denominator needs to be evaluated just once for all pixels and (2) a "map" of the 
numerator is simply an inverse Fourier transform (IFT). Therefore, the entire filtering process can 
be accomplished with ~ ©(A'^pixei log2 A^'pjxei) calculations. On a desktop computer with a 2.5 GHz 
Intel processor and 4.7 GB of RAM, the ~ 200,000 pixel map of Fig. [T]can be filtered according 
to Eq. [121 using the IDL programming language, in under 0.5 seconds. 

Using the standard matched filter is similar to performing a simple least squares minimization 
under the approximation that all pixels in the map have the same uncertainty. If the underlying 
noise is not uniform and this method is used, one must recognize that W ^ C~^ and so the error 
distributions of the Sp cannot be formally calculated as given in Eq. [10] above. 



4. The case of non-uniform coverage 

We show here that, in a map with non-uniform coverage, it is possible to treat the general 
problem of point source identification in a computationally efficient, yet mathematically sound way. 
We find that loosening the assumption of uniform coverage results in an increase in computation 
time that is noticeable but inconsequential in practical terms. As the calculation of D turns out to 
be more complex than the calculation of A^, we tackle N first. 



4.1. The calculation of N 

N, as expressed in Eq. [6l can be calculated through a sequence of FTs and IFTs. We motivate 
this by casting Eq. [6] in a more suggestive manner: 



j27rjkaXk 

Proceeding from the innermost parentheses outward, the steps needed to calculate A'^ are: 



-.^^E^%^(E^S;;;^(^E§^— ]) (-) 



1. Form a new map M, filled with values M(xi) = (i(xi)/e(xi). 

2. The calculation within the innermost parentheses can be recognized as a Fourier Transform 
(see Appendix A). Therefore, form M = FT[M]. 

3. Similarly, the calculation within the next set of parentheses is an IFT. Thus, form P = 

IFT[M/y2]. 

4. Form a new matrix Q, which is filled with values (5(xk) = P(xk)/e(xk). The sum over k in 
Eq. [13] is simply a convolution of Q by the PSF /. 

5. As it is computationally advantageous to carry out this convolution in the Fourier domain, 
calculate / = FT[/] and Q = FT[Q]. 

6. Finally, iV = IFT[/*Q]. 

There are 5 FFTs in this sequence. On the same standard desktop PC described in Section [3l the 
complete computation sketched above takes less than one second for a ~ 200, 000-pixel map. 



4.2. The calculation of D 

The calculation of the denominator, D, is more involved than the numerator as there are 
no simplifications that allow D to reduce to a series of Fourier Transforms. As a result, a full 
calculation of D requires of C(-^pixei) calculations. In the case of the ~ 200, 000-pixel map used 
here, this takes about seven hours on the basic computing platform considered in Section [3j We 
can do better though. Below we give an approximation to D with an accuracy better than 0.1% 
that takes ~ 2 minutes to calculate. 

Our Approximation of D makes use of the band-diagonal nature of the noise correlation matrix 
or, equivalently, the smoothly varying nature of ^^(ka) as discussed in section [2j We start by 
defining a new vector, xj = x^ — xi, and re-writing Eq. [7] as 

n_ 1 v^^,„ ^^^ /(xl-Xp+Xd)/(xl-Xp) 



10 



where Z(xd) is the IFT of l/y2(ka)^ or 

To clearly demonstrate that the sum over / is a convolution, we define the new functions F^ and 
Rd as follows: 

F,(x) = /(x + xd)/(x) (16) 

Rdi^) = , ^\ . V (17) 

e(x + Xd)e(x) 

Thus, Ffi (or R^j is the product of / (or 1/e) and a version of / (or 1/e) that is shifted by the 
vector Xd- Using these definitions, we may express Dp as 

^P = l^E^(^d)5]i^d(xi-Xp)i?d(xi), (18) 

pixel d I 

which is recognizable as a convolution. Finally, using the discrete convolution theorem (see Eq. lASp . 
we find that 

Dp = ^ y. ^(xd) ( y F/(ka)i?,(ka)e~2-^'^--p 1 . (19) 



J— j;Z(xd) ( j;F/(ka)i?,(ka)e~ 



The advantage of casting D in the form of Eq. [19] is the following: Z(xd) is narrowly peaked 
near Xd = because it is the IFT of a smoothly varying function 1/1/^ (ka). Therefore, the sum over 
d converges very rapidly at small Xd and changes very little at large Xd and so high accuracy in the 
calculation of D can still be obtained after truncating the calculation at Nt <C A'^pixei terms. The 
choice of Nt will depend on the properties of .Z'(xd) and the desired accuracy in D. For example, 
for our AzTEC maps we find that if we limit the sum over d to the Nt = A'^pixei/200 ~ 1, 000 terms 
where |Z(xd)[ is largest, D converges to within 0.1% of its final form in useful parts of the map 
(the > 5% coverage region) and requires only 2 minutes to compute. 

The steps we follow to generate our approximate D are then 

1. Compute Z = IFT[l/y^] and find the Nt positions Xd where \Z\ is largest. Then, for each of 
those Xd, 

2. shift / by Xd and multiply by the unshifted / to form F^. Then generate R^ by performing 
the same steps on 1/e; 

3. compute the Fourier Transforms F^ = ¥T[Fd\ and Rd = ¥T\Rd\] 

4. compute the term within parentheses in Eq. [19]as G(xd) = WT{Fd Rd)', 

5. Repeat for the Nt vectors (xd) chosen, and sum the terms Z(xd)G(xd)/A^pixei- 



11 



5. Application of the optimal filter to AzTEC maps 

Here, we will demonstrate how the methods developed above have been applied to AzTEC 
maps, using the field of Fig. [J as an example. Although the examples and justifications presented 
here are based on AzTEC data, we note that many of the trends and techniques identified here are 
typical of most observations and analysis chains. 



5.1. Validation of assumptions using the PSF 



In AzTEC, we have always used an accurate simulation of the PSF, rather than a generic 
form such as a Gaussian, as the template /(x) used in the filter. According to the methodology 
developed in section [21 this is the correct template to use, and should lead to higher accuracy 
of filtered maps. The simulation used for generating the PSF includes effects such as individual 
detector beam shapes, the location and orientation of the field during each observation, as well 
as "cleaning" and filtering steps identical to those used on t he true field. The d etails of how the 



PSF i s generated in the AzTEC data analysis can be found in lDownes et al.l (J2012l ) and lScott et al 



(|2008l ). Fig. EJa) shows the PSF that applies to the AzTEC map of Fig. [H Strictly speaking, the 
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Fig. 2. — (a) The PSF resulting from optics and data "cleaning." The grey scale is stretched to 
emphasize the asymmetric and, at times negative, low-signal values. The dark high-signal regions 
follow a rotationally symmetric profile. The three white contours are at 10%, 50%, and 90% of the 
peak value. The peak value is 0.93% of the inserted source brightness, indicating the amount by 
which a point source is attenuated due to the "cleaning" of atmospheric contamination, (b) The 
rotationally symmetrized version of (a) used as the fit template. The central region has a FWHM 
of ~ 29". 

PSF of Fig. ^a) only applies to a point source at a particular location on the map. Therefore, 
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before adopting assumption 1 of section 12.2^ we constructed several PSFs that apply to different 
locations on a map and verified that th ey are all very si milar (at the few percent level) in terms 



of amplitude, width, and general shape (jScott et al.ll2008l ). The only noticeable difference between 



these PSFs is the orientation of the low-level asymmetric features seen in Fig[2]^a). Therefore, we 
use a rotationally symmetrized version of the PSF, shown in Fig. [2]^b), as the filter template. To 
evaluate the inaccuracy introduced by rotationally symmetrizing the PSF, we have fit the PSF 
of Fig. EJb) to the PSF of Fig. [2|^a), and find that the best fit amplitude differs from the actual 
amplitude of Fig. ^a) by ~ 0.5%, which is small compared to other errors we expect. Therefore, 
we adopt assumption 1 and the rotationally symmetrized PSF of Fig. El^b). 

Assumption 2 of section 12.21 would be difficult to adopt if the field is imaged only once, so 
that the field has an essentially fixed Az-El orientation during the observation and each point on 
the map is scanned over in essentially one direction. However, AzTEC and other mm/sub-mm 
surveys routinely adopt assumption 2 because each of the mapped fields is imaged many times with 
varying orientations and scan directions. This is evident in Fig. [2|^a) where the grey scale has been 
stretched to highlight the low-level asymmetric features of the PSF. We believe that such features 
are due to a small anisotropy in the final aggregate of scan directions from the 45 observation of 
this field. However, the fact that these features are so faint compared to the symmetric central 
parts (refer to the color bar of Fig. [2|) supports our adoption of assumption 2. 



5.2. Estimation of the PSD 

In the AzTEC data analysis, the PSD corresponding to an observation is evaluated using 
Fourier transforms (FTs) of noise realization maps. Our noise realizations are generated using 
a "jackknife" method, where the signs of the detector tim e-streams are switched many times for 



each observation prior to map making (jScott et al.ll2008l ). An important consideration is to use 



long-enough time lags between the random sign switchings in order to preserve long-range noise 
features (low-frequency map noise) due to residual atmospheric contamination and instrumental 
drifts. The ~ 10 s time scales used are longer than the scan turn-around times but short enough 
that astronomical signal will not show up in noise realizations. Next, the noise realization time- 
streams are put through the same cleaning and filtering steps as the actual data, before maps are 
made. For a map such as Fig. [H we generate 100 independent noise realization maps. 

Because noise realizations are used for generating the PSD, it is free from astronomical sig- 
nal and only includes the two dominant noise sources, residual atmospheric contamination and 
instrumental noise. The sub-dominant yet noticeable contribution of confusion noise is left out 
from the PSD for convenience!^ However, the method presented here, of using a PSD to char- 



^ We note that other analyses that make use of the fihered signal map, such as estimating de-boosted source fluxes, 
stacking analyses, and number- counts estimates will not be skewed by the omission of c onfusion noise due to the use 
of source realizations and/or the AzTEC PSF having zero mean fsee lScott et al.ll2008l l. On the other hand, source 
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acterize noise correlation s, can fully accomm odate astronomical "noise sources" such as confu- 
sion noise. For instance, IChapin et alj (|201lh describe how confusion noise was included in their 
matched filter. In general, astronomical effects such as confusion noise are frequency dependent. 
At longer wavelengths, a method fo r including the CMB within a matched filter is described in 
Tegmark fc de Oliveira-Costal (|l998l l. 



Scale [arcsecl 

100 
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Fig. 3. — The average power spectral density (PSD) of noise in noise-realization maps. The broad 
peak at the PSF scale (~ 30") is believed to be caused by atmospheric emission that is imaged in 
between the random sign switches used to generate noise realizations. The dotted line indicates an 
alteration to the PSD (of plateauing at low frequency) that we have tried. 

In our case, the noise PSD is generated by truncating each noise map to include only the 
central > 70% coverage region (see Fig. mb)). Then, the 2-d FTs obtained from these regions are 
rotationally symmetrized and averaged to obtain the PSD shown in Fig. [3l As expected, the noise 
increases with decreasing k over most of the k range. A broad peak is visible at ~ 30", which is 
close to the FWHM of the PSF. This, we believe, is caused by the optical imaging of atmospheric 
fluctuations in between the random sign switches used to generate noise realizations. Near the 
~ 100" scale, the PSD turns over and starts to decrease. This decline cannot be a map-size effect 
because the > 70% coverage region used to construct the PSD extends > 700" in all directions. We 
believe that the PSD's turn-over is real and that it is caused by the principle component analysis 
(PCA) based cleaning that the data is put through prior to map making, in order to mitigate 
the effects of atmospheric contamination and detector cross-talk. As PCA cleaning makes use of 



lists for follow-up observations, which are based purely on s/n will be very slightly biased due to the omission of 
confusion noise. This bias would favor peaks found in regions that have an over-density of unresolved sources, which 
can be beneficial to our understanding of SMGs and their environments. 
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detector-detector time-stream correlations to "subtract" out these effects, it makes sense that the 
PSD shows a decUne on the scale of the detector array (the spacing of individual detector beams is 
~ 40" and the footprint of the entire array is ~ 480" for AzTEC-ASTE). To test if this long-range 
decline in the PSD is indeed real, we have performed the optimal filter presented here using the 
measured PSD as well as a modified PSD that plateaus at small k, as indicated by the dotted line 
of Fig. [31 The fact that there is no perceptible difference in results indicates that the map noise 
power is indeed low at low k. As noted earlier, an important property of the PSD of Fig. [3] is that 
it is a smoothly varying function, as opposed to a narrowly peaked one. 



5.3. Filtering of maps 

Once the PSF and PSD are available, the methods of section [3] can be applied to obtain the 
filtered signal and coverage maps. In practice, more time is spent on the calculation of N than the 
calculation of D during the filtering process, even though D involves a more complex calculation 
(see section U]). This is because N needs to be evaluated separately for the actual map and each 
noise realization, as N depends on d{x.) (see Eq. [6]), while D needs to be evaluated only once for all 
101 maps. For the example field used, the filtered signal and coverage maps are shown in Fig. [H 





1b7.8 1b7.7 157.6 157.5 157.' 
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(b) 



Fig. 4. — Filtered signal map and filtered coverage map of the entire field, generated according to 
the optimal filter of section HI 



The filtered coverage map must be filled with values I /rip, and therefore, according to Eq.[TOl it 
is simply equal to the map D. But D is generated by propagating the initial noise variances, e(xp), 
according to Eq. [71 Thus the noise estimates contained in D are only as good as the starting noise 
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estimates e(xp), which are computed from estimates of individual de tectors' noise leve ls "near the 



time" that they contribute a sample to the map pixel in question (see lScott et al.ll2008l . for details). 
Therefore, before proceeding to the step of identifying point sources, we check the accuracy of noise 
estimates contained in D using noise realization maps. First, we form a second noise estimate n' 
for each pixel by taking the standard deviation of the 100 noise-map values found at that pixel. 
Next, we form a second coverage map T filled with the values 1/n'p . Although T provides a robust 
estimate of coverage, it is noisier than the the original coverage map D due to the sample size of 
100 (we expect a ~7% error in Tp for this sample size). Despite the noisy appearance of T, its 
overall shape agrees very well with D. However, T and D often differ from each other by an overall 
scaling factor. We determine this scaling factor by comparing the average values of T and D in the 
>70% coverage region and applying a correction factor to D in order to obtain the final filtered 
coverage map. Thus, the filtered coverage map is rjD, where 

V=^, (20) 

and 7] usually lies in the range 0.85-0.95, depending on the field. Although we use just the >70% 
coverage region to find rj, the agreement between T and rjD remains good out to the very edge of 
the mapped region until the coverage dips below 5%. Beyond this region, rjD is consistently larger 
than T. Therefore, we do not trust the filtered coverage map beyond the 5% region and exclude 
this area when searching for point sources. Finally, in Fig. [5l we present the filtered signal-to-noise 
map, which is used directly to find point sources. It is generated by dividing the filtered signal map 
(Np/Dp) by the filtered noise map (l/y^r/Z)p)o 



6. Results and discussion 

It is encouraging that the signal-to-noise map of Fig. [5] has its most prominent peaks near 
the center, where the coverage is highest, and that the number of peaks as well as their amplitude 
declines smoothly toward the outer regions. Thus, even though the filtered signal map has large 
fluctuations near the edges (see Fig.l^a)), the noise in those regions are accounted for in this optimal 
filter and there in no need for a by-hand coverage cut to eliminate spurious behavior near the edges. 
Throughout the AzTEC data analysis campaign, which has led to the publication of many source 
lists, variations of the filter presented here have been used. Therefore, the effectiveness of this 
technique may be assessed by the success rate of follow-up observations. In this regard, AzTEC 
has enjoyed an excellent record thus far. For instance, of 15 Az TEC point sources fo l lowed up by 



the Submillimeter Array (SMA), all were successfully observed ([Younger et al.l 120071 . l2009l ). This 



is a very good record in comparison to the usual success rate of mm/sub-mm follow up studies. 



^The correction factor r] does not need to be applied in generating the filtered signal map because A'^ and D both 
contain two factors of e (see Eq. |6]and Eq. [7]). 



16 




57.8 157.7 157.6 

RA [deg] 



157.5 



157.4 



s/n 



Fig. 5. — The filtered signal-to-noise map. The 5%, 30%, and 70% coverage contours are overlayed. 



In terms of weighting a fit according to all possible sources of uncertainty, the filtering technique 
presented in section H] (filter A) is clearly closer to optimal than the standard matched filter discussed 
in section [3] (filter B), especially in the presence of significant coverage gradients. However, filter A 
contains a larger number of steps and computations. In terms of implementation, we have shown 
that the intricacies of filter A have virtually no effect on the total computing time devoted to a 
field. However, converting these additional steps into functional computer code does require more 
time and effort. Therefore, it is fair to ask how much better filter A performs compared to filter B 
in a practical application. Of course, the answer depends on the specifics of the application. For 
instance, the acuteness of coverage gradients in a map will be important for this assessment. We 
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will try to answer this question in the context of the typical AzTEC map that we have used here 
for illustrative purposes. 

First, when a fit is closer to optimal, it yields smaller error bars on the fit parameters than a 
less optimized fit. A fair way of comparing error bars here is to compare the T maps (generated 
using noise realizations, as described in section [53]) that result from the two filters. We find that 
the post-filter coverage (l/Up) is on average 5-6% better with filter A than with filter B over most 
of the map, including the > 70% coverage region. Thus, filter B is more prone to error, which in 
this context leads to higher rates of false detections as well as non-detections of truly significant 
sources. 

To form an impression regarding these rates, we summarize in table [U the results obtained with 
each filter in different regions of the map. The "source candidates" reported in table [1] are defined 
to be peaks in the signal-to-noise map (Fig. [5]) that exceed a value of 3.5. This is a good choice of 
threshold because, on average, pure noise realizations yield only ~ 1 noise peak with signal/noise 
> 3.5 over the entire map. Given the non-uniformity of coverage within the 30-70% region, it is 
not surprising that filter B "detects" 3 peaks that are judged by filter A (more accurately) to have 
s/n < 3.5. It is interesting that even in the > 70% coverage region, the source lists generated with 
filters A and B differ by two sources. The results reported in table [T] provide some grounds for the 
reader to gauge whether the additional work involved in implementing filter A is worth the effort. 
Of course, the optimal filter presented here will be most useful when the coverage of a field has 
large non- uniformities, unlike in the example field used here. 

Thus far, many AzTEC publications only provide source lists from map regions with > 70% 
coverage. As the optimal filter presented here is geared to handle coverage non- uniformities, we 
believe that the reliable region for source extraction can be extended to lower coverage thresholds. 
As mentioned in section 15.3] regions with < 5% coverage should be excluded. In fact, it may be 
best to exclude a larger region. For instance, not all elements of the detector array have imaged 
the outer regions of the map. Therefore, those regions may be systematically biased in some way 
compared to regions viewed and appraised by all detectors. Given that the AzTEC- ASTE detector 

Table 1. The number of peaks with s/n > 3.5 found with the fully optimal filter of section [H 

(filter A) and the standard matched filter of section [3] (filter B) are reported. The last column 

gives the number of "sources" found with filter B but not with filter A. The 30-70% coverage 

region has an area that is 75% of the >70% coverage region. 



Map region Filter A Filter B B but not A 

>70% coverage 31 31 2 

30-70% coverage 15 18 3 
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array has a foot print on the order of 480", it is reasonable to leave out a border of half that size, 
and use the > 30% coverage region of Fig. [5] for source identification. This represents an increase 
of 75% in the source extraction area. Table [1] shows that ^ 50% more source candidates can be 
obtained this way. 

The optimal filter described here is now part of the AzTEC data analysis pipeline, which may 
be downloaded (by anonymous ftp) 



at http://www.astro.uniass.edu/aztec/Software/software.html Within this suite, the par- 



ticular IDL routine that implements the optimal filter 

is aztec-adapative-wiener_filter.pro. 

This work has been funded, in part, by NSF grant AST-0907952. KSS is supported by the 
National Radio Astronomy Observatory, which is a facility of the National Science Foundation 
operated under cooperative agreement by Associated Universities, Inc. 



A. Appendix: Fourier transform conventions 

In this work, FT[(7] denotes the discrete 2-dimensional Fourier transform of the function (?(x). 
Following the convention used by IDL and other high-level programing languages, this Fourier 
transform is defined as 

ff(k) = gikj + hj) = — — Y. Y. 9ixn.i + y„j)e-2-i«WA^^e-2-^■^"/^^ (Al) 

" m=0 n=0 

where N^ and Ny are the number of pixels in the x and y dimensions. Therefore, N^Ny = A^pixei- 
The Xm and yn above, take on the x and y values of all the pixel centers. The a and h above 
assume integer values in the range [0, Nx — 1] and [0, Ny — 1] respectively. The discrete points on 
the reciprocal plane, where g is defined are specified by 

where Ax = Ay is the pixel size. The corresponding inverse Fourier transform, denoted IFT[^], is 
defined as 

N:,-lNy-l 

g{^)=g{x^l + ynj) = Y H ~g{kai + hj)e^^^'^^"'- e^^^^'^^y . (A4) 

a=0 fe=0 

According to these conventions, the discrete convolution theorem takes on the form 

Y H^i - xp)5(xi) = iVpixci Y /i*(ka)5(ka)e2"^■'^-"^ (A5) 
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when h and g are real functions. In Eq. IA51 / is an index over all pixels of the map rather than a 
single dimension and a, similarly, is an index over all sampled points on the 2-d reciprocal plane. 
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